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At sufficiently low temperatures, condensed-matter systems tend to develop order. An ex- 
ception are quantum spin-liquids, where fluctuations prevent a transition to an ordered state 
down to the lowest temperatures. While such states are possibly realized in two-dimensional 
organic compounds, they have remained elusive in experimentally relevant microscopic two- 
dimensional models. Here, we show by means of large-scale quantum Monte Carlo simu- 
lations of correlated fermions on the honeycomb lattice, a structure realized in graphene, 
that a quantum spin-liquid emerges between the state described by massless Dirac fermions 
and an antiferromagnetically ordered Mott insulator. This unexpected quantum-disordered 
state is found to be a short-range resonating valence bond liquid, akin to the one proposed 
for high temperature superconductors. Therefore, the possibility of unconventional super- 
conductivity through doping arises. We foresee its realization with ultra-cold atoms or with 
honeycomb lattices made with group IV elements. 

The quantum mechanical description of the relativistic electron was attained by Dirac, who 
revealed both its intrinsic angular momentum (the spin), with a half-integer quantum number 
S = 1/2, and the existence of its antiparticle, the positron'. Both obey the Fermi-Dirac statis- 
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tics, which implies that two identical particles cannot occupy the same quantum mechanical state. 
Such particles are generically called fermions. In case of a vanishing rest mass, the energy of 
Dirac fermions is a linear function of momentum. Such massless Dirac fermions were recently 
observed in two-dimensional solid-state systems like grapheneP^^ and surfaces of bismuth based 
compoundd^. For graphene, a single layer of carbon atoms with honeycomb structure, unusual 
electronic behaviour is anticipated, and partly verified experimentally, due to the two-dimensional 
Dirac-like dispersion of the electrons at low energieP. The interplay of a relativistic dispersion 
with interactions at half-filling is expected to lead to a quantum phase transition between the 
semimetal (SM) at low and a Mott insulator (MI) at high interaction strengthd^. Here, a Mott 
insulator is an insulating state that results not from the band structure alone, but is due to the ef- 
fects of interactions. Such correlation effects can be displayed by the Hubbard model in its most 
basic form, as exemplified in high temperature superconductivity^, or with ultra-cold fermionic 
atoms loaded in optical lattice^^'^] Studies of Hubbard-like models on the honeycomb lattice 
suggested the emergence of exotic phases such as gapless spin liquid^'^, charge density wave^, 
quantum spin Hall stated^^CS^ or superconductivity^^ at or near a density of one fermion per site 
(half- filling for the two-species case). 

Given the various phases proposed for fermions on a honeycomb lattice based on Hubbard- 
like models, it is important to explore the ground-state properties in the intermediate coupling 
regime of the original lattice model with an unbiased method. Due to the absence of a sign-problem 
in determinantal quantum Monte Carlo (QMC) simulations (see Methods) in the half-filled case, it 
is the method of choice for extrapolations to the thermodynamic limit (TDL), leading to essentially 
exact results limited only by the statistical noise. Employing large-scale quantum Monte Carlo 
simulations of the spin-| Hubbard model at half-filling on the honeycomb lattice, we show that for 
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intermediate interactions a gapped non-magnetic phase destroys the semimetal before the transition 
to an antiferromagnetically ordered Mott insulator at strong interactions sets in. This quantum 
spin-liquid phase is characterized by local correlations that correspond to a resonating valence- 
bond (RVB) stateP^Elas proposed in the context of high temperature s uperconductivit)^^*^^*^. 

Following their original proposal^^^'t^, spin-liquid states were established in effective mod- 
els of singlet-dynamics such as quantum dimer modelP^^. Our results show that RVB states are 
realized in a microscopic model of correlated electrons, bringing closer their observation in ex- 
periments. Honeycomb lattices of group IV element^ and ultra-cold fermionic atoms loaded in 
optical lattice^'^'^ appear as promising candidate systems to realize the RVB state out of Dirac 
fermions. 

Phase diagram from quantum Monte Carlo 

Previous numerical studies of the Hubbard model on the honeycomb latticeP^I^ suggested that 
a single quantum phase transition separates the paramagnetic weak-coupling SM phase from a 
strong-coupling antiferromagnetic (AF) MI. At strong enough repulsion, antiferromagnetism is 
certainly possible since the honeycomb lattice is bipartite, so that AF order is not geometrically 
frustrated. However, the honeycomb lattice has the smallest coordination number in two dimen- 
sions, such that the effect of quantum fluctuations is the strongest. Hence, the competition between 
the tendency to order and quantum fluctuations requires a detailed analysis of correlations and a 
careful extrapolation to the TDL in order to characterize the possible phases. Here, we present 
results based on projective (temperature T = 0) determinantal QMC simulations in the canonical 
ensemble at half-filling. In order to assess the above scenarios, we focus in particular on the region 
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near the Mott transition. 



The Hamiltonian of the spin-| 



i Hubbard model on the honeycomb lattice equals 
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where c-^ (Cja) denotes the creation (annihilation) operator for fermions of spin a =t, i on lattice 
site i, and riia = c\^Cia. Here, t denotes the nearest-neighbour hopping amplitude, and > the 
strength of the onsite repulsion. Our notations in real and momentum space are shown in the inset 
of Fig. 1. At f/ = 0, the tight-binding Hamiltonian has a linear dispersion near the Dirac points 
{K, K' - cf. Fig. 1), where the conduction and valence bands touch at half-filling, corresponding 
to a density {uia) = 1. At half-filling, the finite-t/ region can be studied using projective QMC 
to obtain ground-state expectation values of any physical observable. Details are relegated to the 
Methods section. The phases described in the following correspond to extrapolations to the TDL. 
For that purpose we study lattices of = sites with periodic boundary conditions, and linear 
sizes up to L = 18. 

To monitor the electronic properties of the system upon increasing U, we extracted the single- 
particle excitation gap Asp(k) from the imaginary-time displaced Green's function (cf. Supple- 
mentary Information (SI) for details). Asp(k) gives the minimal energy necessary to extract one 
fermion from the system, and corresponds to the gap that can be observed in photoemission ex- 
periments. As shown in Fig. 1, Asp(A') = for t/ /t below about 3.5, as expected for a SM. For 
larger U /t, the system enters into an insulating phase due to interactions. The values of the gap are 
obtained via an extrapolation of the QMC data to the TDL as shown in Fig. 2a. 

From previous analysis of the model, one expects long-range AF correlations when the MI 
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appears. We therefore measured the AF spin structure factor Saf (cf. SI) that reveals long-range 
AF order if = limAr„j.oo Saf/N > 0. Figure 2b shows the QMC results together with a finite 
size extrapolation. The results of the latter are also presented in the phase diagram of Fig. 1. AF 
order appears beyond = 4.3, a value that is consistent with previous estimates for the onset of 
long-ranged AF ovde^^. This leaves an extended window 3.5 < U/t < 4.3, within which the 
system is neither a SM, nor an AF MI. 

Further details on the nature of this intermediate region are obtained by examining the spin 
excitation gap, extracted from the long-time behaviour of the imaginary-time displaced spin-spin 
correlation function (cf. SI). We consider first the spin gap in the staggered sector at k = F, 
which vanishes inside the AF phase due to the emergence of two Goldstone modes, as well as in 
the gapless SM phase. Figure 2c shows finite size estimates of for different values of U/t, along 
with an extrapolation to the TDL. A finite value of A^ persists within an intermediate parameter 
regime 3.5 < U/t < 4.3, while it vanishes both within the SM and the AF phase. This dome 
in the spin gap is also seen in the inset of Fig. 2c, that displays both the finite-size data and the 
extrapolated values of A^ as a function of U/t. We also calculated the uniform spin gap A^, by 
extrapolating the spin gap observed at the smallest finite k-vector on each cluster to the TDL. 
Au is found to be even larger than A., inside the intermediate region (e.g. A„ = 0.099 ± 0.001 
(s.e.m.) at f//t = 4), and vanishes in the SM and the AF phase (A„ cannot be measured directly at 
k = 0, because the uniform magnetization is a conserved quantity, cf. SI). Hence, this intermediate 
insulating region corresponds to a spin-gap phase. 

From analysing the [/-dependence of the kinetic energy density, iJfcjn = (— t^^j^-^ 
Cj^Cia))/N, we obtain further insight into these different regimes and the emergence of local mo- 
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ments. As shown in Fig. 3, the curvature (fEkin/dU'^ changes sign near U/t = 4.3. This marks 
a characteristic change from the weak-coupling region of positive curvature with delocalized elec- 
trons to the strong-coupling AF region with negative curvature. In the latter region, localized spins 
form and order in an AF state. In the intermediate spin-gap region, fluctuations are large enough 
to still prevent the formation of well-localized magnetic moments. Note, that around U/t = 3.5, 
a change in the curvature can be observed, that adds to the already presented evidence for an 
intermediate phase. 

Characterization of the spin-gap phase 

The observation of a finite spin gap rules out a gapless spin-liquicP^^Il, quantum spin Hall stated, 
as well as triplet superconductivit)^^: The remaining possibilities can be enumerated by consider- 
ing the coupling to order parameters that lead to the opening of a mass gap in Dirac fermions-^, and 
hence account for the single-particle gap observed in the QMC data: (i) singlet superconductivity, 
(ii) a quantum Hall state (QHS j^, (iii) charge density wave (CDW) ordei^ and (iv) a valence bond 
crystal (VBC). 

In order to assess if superconductivity arises in the vicinity of the Mott transition, we used the 
method of flux quantisation which probes the superfluid density and is hence independent of the 
specific symmetry of the pair wave functiorP^I. Let $ be a magnetic flux traversing the centre of a 
torus on which the electronic system lies and Eo{^ / <l>o) the total ground-state energy, $o being the 
flux quantum. A superconducting state of Cooper pairs is present if in the TDL, the macroscopic 
energy difference -E'o($/$o) — Eo{^/^o = 1/2) is a function with period 1/2!^. In contrast, a 
metallic phase is characterized by i?o($/$o) — Eo{(^/^o = 1/2) vanishing as a power law as a 
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function of system size, while in an insulating phase, it would vanish exponentially. As shown in 
the SI, this quantity vanishes in the TDL both in the semi-metallic state at t/ = and at = 4, 
i.e. in the intermediate phase. In addition, we measured pair correlations, ruling out superconduc- 
tivity in (extended) s-, p-, d-, and /-wave channels (cf. SI). Hence, both flux quantization as well 
as a direct measurement of pair correlations lead to no sign of superconductivity. 

Both the CDW and QHS trigger a breaking of the sub-lattice symmetry and thereby open a 
mass gap at the mean-field level. A detailed analysis of the charge-charge correlation functions 
rules out a CDW. Furthermore, we find no signature for the presence of (spin) currents in the 
ground-state (cf. SI). This rules out the breaking of sublattice and time reversal symmetries, as 
required for the QHS, in the pristine Hubbard model, and possibly, extensions of it are necessary 
to reach such a stateP'. 

To examine the occurrence of a VBC, we probed for dimer-dimer correlations between sep- 
arated dimers formed by nearest neighbour bonds (ij) and (kl) (cf. SI). We find no VBC, neither 
in the charge, nor in the spin sector. Figure 4 shows the results of this measurement in the spin 
sector, i.e. the correlation between singlet dimers at U/t = 4.0. The striped bond is the one with 
respect to which correlations were determined. They are found to be short-ranged, and consistent 
with the dominance of a RVB state within the hexagons of the honeycomb lattice. This can be seen 
by comparing the singlet-correlations with those of an isolated hexagon (inset Fig. 4), the classi- 
cal example of the resonance phenomenon in conjugated 7r-electrond^. Accordingly, we find no 
long-ranged order from the dimer-dimer structure factors in Fourier space. Our results thus reveal 
a genuinely exotic state of matter, where no spontaneous symmetry breaking is observed, while a 
spin-gap is present. It corresponds to a spin-liquid RVB state in the intermediate coupling regime 
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in the vicinity of the Mott transition. 
Further insight into the RVB state 

The QMC results presented above uncover the realization of a quantum spin-liquid state of corre- 
lated fermions on a non-frustrated, bipartite lattice. In principle, such quantum-disordered states 
can occur in different flavors, and we thus aim to shed further light on the nature of the spin-liquid 
observed above. Gapless (algebraic) spin-liquids, or long-range RVB states are characterized by 
critical spin-spin correlationd^llMlSl observation of a finite spin-gap clearly excludes such can- 
didate states, while being consistent with the characterization of the observed spin-liquid in terms 
of short-range RVB state^El'^^'^, also in accordance with the observed short-ranged dimer-dimer 
correlations. Short-range RVB states are modeled in general by quantum dimer modeld^^^, which 
capture the fluctuations of singlets in a RVB state, with dimers being a strong-coupling represen- 
tation of nearest-neighbour singletd^. Depending on the lattice geometry, quantum dimer models 
can exhibit a fully gapped short-range RVB phase, as e.g. on the triangular lattice!^, but also ex- 
hibit spin liquid states with gapless excitations, as is the case for bipartite lattices at the Rokhsar- 
Kivelson point^' and within the U(l) spin-liquid phase stabilized for spatial dimensions d > 
In the later case, gapless singlet excitations constitute an emerging "photon" soft mode. Fully 
gapped phases of quantum dimer models are furthermore characterized by a non-trivial topolog- 
ical order, implying e.g. an emerging ground-state degeneracy of two-dimensional systems with 
periodic boundary conditions in the TDL'^'^. 

In order to assess, whether topological order can characterize the short-range RVB in our 
case, we examine the low-energy singlet excitations. As proven in an exact theorem by LiebPl, 
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the finite systems used in our numerical simulations have a non-degenerate singlet ground-state for 
any finite value ofU > 0. Hence, degeneracy can only appear in the TDL. If so, low-energy singlet 
states should be present, with decreasing excitation energy as the system size increases. Since our 
QMC method projects out the finite system's ground-state from a singlet trial wavefunction, we 
can monitor the expectation value of the internal energy E{Q), where 6 is the projection parameter 
(cf. Methods). Here all contributions from singlet states with the same quantum numbers as the 
ground-state are included, that have a finite overlap with the trial wave function. For a given system 
of size N, we define 9* as the value of the projection parameter, such that {E(Q) — Eq) /N < e, 
for 9 > 9*, where Eq is the ground-state energy, and e is an energy scale which we choose small 
enough, in order to guarantee that e <^ {Ei — Eq)/N, where Ei is the lowest singlet excited 
state above the ground-state. Typically, we choose e of the order of our statistical error in the 
energy density. Then, 1/9* is a lower bound for the lowest singlet excitation (cf. SI). Such an 
analysis on various lattice sizes leads to the conclusion, that the singlet excitation gap stays above 
the spin gap in the TDL (cf. SI), providing no evidence for the emergence of a topological state. 
However, we cannot definitely exclude topological order, if the relevant singlet states happen to 
have a vanishing overlap with our trial wave function. For the future, it will be interesting to 
explore the low-energy singlets beyond the projective scheme, and probe for soft modes similarly 
as does e.g. the construction of finite momentum trial states in quantum dimer models 1^^^. Our 
findings, based on a controlled numerical framework, therefore open a new facet of quantum spin- 
liquids, where an appreciable amount of doubly-occupied sites are present, extending well beyond 
the regime of localized spin physics. 
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Discussion and outlook 

The presence of a spin-liquid in the Hubbard model on the bipartite honeycomb lattice close 
to an antiferromagnetic Mott insulator resembles the situation in the organic antiferromagnet 
fi:-(BEDT-TTF)2Cu2(CN)3, which has been argued to display a spin-liquid stateP^ED, albeit the lat- 
ter system is on a triangular lattice and hence frustrated. This difference can be reconciled starting 
from the strong-coupling limit of the Hubbard model, i.e. a nearest-neighbour Heisenberg model, 
that close to the Mott transition acquires corrections that induce efficient frustrations to the spin 
degrees of freedom. In fact, a Klein Hamiltonian for a spin-liquid state on the honeycomb lattice 
was constructed, including extended exchange interaction^^. A more pronounced difference is 
the appearance of superconductivity in the organic systems upon pressure, that is equivalent to a 
reduction of the ratio U /t'm the Hubbard modeP^. The absence of superconductivity in our system 
could be due to the vanishing density of states at the Fermi energy. In this case, a finite coupling 
strength is needed, at least in the BCS-frameP'. However, having an unexpected realisation of a 
short-range RVB state, it would be highly interesting to explore the consequences of doping, in 
a spirit rather close to the original scenario proposed by AndersonP' and Kivelson et al^ for the 
cuprates. In particular, for the fully gapped short-range RVB state, the finite spin-gap sets the 
energy scale of pairing in the superconducting stateP^. In this respect, the value obtained for the 
spin-gap is rather promising. The largest value attained is ~ 0.025t (Fig. 1), that for t in the 
range of 1.5 to 2.5 eV (in graphene is t = 2.8 eV% corresponds to a temperature scale ranging 
from 400 to 700 K. 

Although studies of doping are beyond the power of our quantum Monte Carlo approach due 
to the sign problem, they could open interesting perspectives e.g. in future experiments with ultra- 
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cold atoms on a honeycomb optical lattice, or with honeycomb lattices based on group IV elements 
like expanded graphene (to enhance the ratio U/t) or Si, where the nearest neighbour distance is 
expected to be approximately 50% larger than in graphenepl, such that correlations effects are 
enhanced. In fact, first attempts succeeded in synthesizing single-crystal silicon monolayers!^. 

Methods Summary 

At half-filling, the finite-t/ region can be studied using the determinantal projective QMC algo- 
rithm to obtain ground-state expectation values of a physical observable by performing an imagi- 
nary time evolution of a trial wave function that is required to be nonorthogonal to the ground-state. 
The value reached in the imaginary time evolution corresponds to a projection parametei^^^. 
For a spin-singlet trial wave function, we found = 40/t to be sufficient to obtain converged 
ground-state quantities within statistical uncertainty. In the presented simulations, we used a finite 
imaginary time step Ar = 0.05/t. We verified by extrapolating Ar — )■ that this finite imag- 
inary time step produces no artifacts. Simulations were performed for systems of linear size 
L = 3, 6, 9, 12, 15 and 18 with = sites. For periodic boundary conditions these clus- 
ters all have nodal A' -points and hence allow a smooth extrapolation to the TDL. Imaginary time 
displaced quantities are obtained by using the approach in Ref. 50. 
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Figure 1 | Phase diagram for the Hubbard model on the honeycomb lattice at half- 
filling. The semimetal (SM) and the antiferromagnetic Mott insulator (AFMI) are sepa- 
rated by a gapped spin liquid (SL) phase in an intermediate coupling regime. Asp{K) 
denotes the single-particle gap and the spin gap. denotes the staggered mag- 
netization whose saturation value is 1/2. Error bars indicate the standard error (s.e.m.). 
The inset shows the honeycomb lattice with primitive vectors ai, aa, and reciprocal lattice 
vectors bi, b2. Open (full) circles for sublattice A [B), the Dirac points K and K', and the 
M and r point are indicated. 
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ML 



Figure 2 | Finite size extrapolations of the excitation gaps and the antiferromagnetic 
structure factor, a, Single-particle gap at the Dirac point Asp{K) for different values of 
U/t, linear in 1/L. Asp{K) is obtained by fitting the tail of the Green's function (inset) 
to the form e~''^"'^^\ b, Antiferromagnetic structure factor Saf for various values of U/t 
using 3rd order polynomials in 1/L. AF order appears beyond U/t = 4.3, as seen in the 
histograms from a IVIonte Carlo bootstrapping analysis (inset), c, Spin gap at different 

values of U/t, using 2nd order polynomials in 1/L. Error bars in a, b, and c indicate the 
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standard error (s.e.m.). 




u/t 

Figure 3 | Derivative dEkin/dU of the Itinetic energy density as a function of U/t 
for systems of different sizes. The dashed line is a fit to the low-f/ behaviour. The 
inset shows the QIVIC data for the l<inetic energy density Ekm from which the derivative 
is obtained by numerical differentiation. Statistical errors (s.e.m.) are smaller than the 
symbol size. 
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-0.0267 



0.0267 



Figure 4 | Real space plot of the spin dimer-dimer correlations. The correlation func- 
tion Dij^ki (cf. SI) for a L = 6 system at f//t = 4, together with the same correlation for the 
isolated Hubbard hexagon also at f//t = 4 (inset). The reference bonds are dressed with 
stripes. Numbers in parenthesis indicate the standard error (s.e.m.) of the last digit. 
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Methods 



The projective QMC algorithm employed for the simulations presented in this article constitutes 
an unbiased, controlled and numerically exact method which is described in detail in Refs. 47,49. 
Within this scheme, ground-state expectation values of a physical observable A are obtained from 
performing an imaginary time evolution 



Here, we used the fact that the ground-state of the Hubbard model on the honeycomb lattice is 
non-degenerate on any finite lattice at half-filling and for an equal number of sites within the two 
sublatticed^, and implicitly assumed that the trial wave function, |\1/t), has a finite overlap with 
this ground-state. The standard implementation of the algorithm requires the trial wave function to 
be a single Slater determinant. 

The efficiency of the projective approach strongly depends on the choice of the trial wave 
function. To generate optimal trial wave functions, different approaches can be employed. One 
possible strategy consists in optimizing the overlap of the trial wave function with the finite sys- 
tem's ground- state"^*^. Alternatively, one can specify a series of good quantum numbers that char- 
acterise the ground-state. The trial wave function is then constructed as to share the same quantum 
numbers. We have chosen the latter approach for our simulations, and generated the trial wave 
function from the non-interacting tight binding model on a torus, through which we thrust a mag- 
netic flux $, corresponding to the vector potential A = <l>^. In particular, we employed a trial 
wave function of the form |\E't) = |\l/T)t ® \'^t)i, where \'^T)a denotes the ground-state of the 
single particle Hamiltonian in the spin-flavor a Hilbert subspace. 



{A)= lim (^r|e-®^/'Ae-®^/2|^r)/(^T|e-®^|*r). 
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Here, $o = /^e/c denotes the flux quantum. At $ = and for the considered finite lattices of linear 
size L = 3n (n E N), the half-filled ground-state wave function of the above Hamiltonian is 
degenerate. Imposing an infinitesimal twist (we verified that taking $/$o = 0.0001 is sufficiently 
weak) lifts the two-fold degeneracy of the single particle states at the Dirac points K and K'. The 
thereby produced filled shell configuration guarantees the absence of a negative sign problem, and 
has total momentum K + K' = G (G being a reciprocal lattice vector) and zero total spin. We 
used this trial wave function for our zero-flux quantum Monte Carlo simulations. Employing this 
trial wave function, we found an imaginary time projection parameter 6 = 40/t to be sufficient to 
obtain converged ground-state quantities within the statistical uncertainty. 

For the presented simulations, we used an SU(2) symmetric, discrete Hubbard-Stratonovich 
transformation which allows for a direct generalization of the simulation scheme to SU(N) sym- 
metric models^^. In this approach, after performing the standard Trotter- Suzuki decomposition!^, 
the interaction part of the imaginary time evolution operator is expressed as 

l=±l,±2 

with the two functions 7 and t] of the four- valued auxiliary field / = ±1, ±2 taking on the values 

7(±1) = 1 + 76/3 , r]{±l) = ±-^2(3-76) , 
7(±2) = 1 - 76/3 , r/(±2) = ±^2 (3 + 76) . 

The advantage of this representation is the fact, that for each Hubbard-Stratonovich configuration, 
the SU(2) spin symmetry of the Hubbard model is conserved explicitly. The above Hubbard- 
Stratonovich transformation produces an overall systematic error proportional to Ar^ in the Monte 
Carlo estimate of observables which, in comparison to the Trotter error of order Ar^, is however 
negligible. We employed a finite imaginary time step At = 0.05/t and verified upon extrapolating 
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At — 0, that this value produces no artifacts. In order to extract the gaps to the various excitations 
of the system, we calculated in addition to equal-time correlations also imaginary time displaced 
correlation functions. To efficiently calculate these imaginary time displaced quantities, we used an 
approach that was introduced in Ref. 50, and which accounts for the fact, that for a given auxiliary 
field configuration the equal-time Green-function matrix is a projector. 

Finally, we have confirmed the validity of our implementation against exact diagonalization 
results on both L = 2 and L = 3 lattices. 
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Supplementary Information 

In these supplementary sections, we provide further details about the simulation results mentioned 
in the main text. As a convenient notation, in the following cj^.^^ and cl^^^ (CxAa and CxBa) denote 
creation (annihilation) operators for fermions of spin a =t or I, on the lattice site that belongs 
to the sublattice A and B respectively, within the unit cell at position x. Furthermore, nxaa = 
4aaCxaa ^nd T^xa = '^xaa dcnotc the local dcnslty operators, and Sxa = |4aa'^a/3Cxa/3 the 
local spin operators, where cr = (a^, (Jy, a^) is the vector of Pauli matrices and a E {A, B}. The 
corresponding operators in momentum space are obtained from 



where Xyi = (0, 0) and = (0, a), with a the distance between neighbouring lattice sites. Sim- 
ilarly, Fourier components rikaa> '^ka and Ska of the density and spin operators are defined. For 
the following, it is also convenient to introduce the three lattice vectors related to the three nearest 
neighbour bonds, 

5i = (0,0), ^2 = -a2, 53 = ai-a2, 
where ai and a2 are shown in Fig. 1, as well as 

^1 — a2, r2 — a2 — ai, — — ai, r4 — — ri, rs — — r2, rg — —r^ 

connecting a given lattice site to its six next-nearest neighouring lattice sites. For the correlation 
between two local operators Oi and O2, we employ a short notation for the cummulant. 



Most of the following results concern the intermediate spin liquid phase, and we present in those 
cases quantum Monte Carlo data for the representative value of = 4. 




(2) 



X 



((O1O2)) := (Oi02)-(Oi)(02). 



(3) 
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1 Green's function and single-particle gap 

To probe the single-particle properties, we measured the imaginary-time displaced Green's func- 
tion 

G(k, r) = ^Yl (4at(^)Ckat(0)) = IY1 i^lalirWam), (4) 
a a 

where c^l^i'^) = e''^c[Jj^e~'^^. The single-particle gap Asp(k) is obtained from G{k,r) oc 
exp(— rAsp(k)) at large imaginary time r, and corresponds to the particle (or hole) excitation 
energy with respect to the chemical potential /i = at half-filling in this particle-hole symmetric 
system. At U = 0, the single-particle gap vanishes at the Dirac points K and K' (cf. Fig. 1 for 
our notation in momentum space), and we thus considered Asp{K) in detail. The quantum Monte 
Carlo data for G{K, r) and Asp{K) is presented in the main text. 

2 Spin correlations and Saf 

The antiferromagnetic order at large values of U/t resides within the unit cell of the honeycomb 
lattice. Hence, the spin structure factor for antiferromagnetic order relates to the staggered spin 
correlations at the T point (cf. Fig. 1 for our notation in momentum space), 

'5AF = (E(SxA-SxB)]ViV)- (5) 

X 

In addition to the above structure factor, we also probed directly the spin-spin correlation functions 

Q^(x,y) = ((Sx,-Sy,)) (6) 

at the largest available distance = ([L/2 + 1] — 1) ai + ([L/2 + 1] — 1) a2 for different system 
sizes, and performed a finite size scaling of both C^'^ = C^'^(0, d^) and Cf'^ = -Cf'^{0, d^). 
A comparison of the scaling of these quantities to Saf is shown for both U/t = A and U/t = 4.5 
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in Suppl. Fig.[TJ exhibiting the consistency between these different approaches to quantify the spin 
correlations in the ground state. 

3 Spin excitation gaps 

The gaps for spin excitations at momentum vector k are obtained from the imaginary-time dis- 
placed spin-spin correlation functions for both the staggered sector, 

Ssik, t) = {{{SUr) - SkB(r)) • (SkA(O) - Ski?(0)))), (7) 

as well as the uniform sector, 

Su{k, t) = (((SkA(r) + Ski?(r)) ■ (Sk^W + SkB(O)))), (8) 

where Sk,a(T) = e^^Sk.aC""^^. Similarly as for the single-particle gap, the spin excitation gaps 
are obtained from ^^(k, r) oc exp(— rA5(k)), and ^^(k, r) oc exp(— rAM(k)) at large imaginary 
time r. The staggered spin gap = Aj,(r) can be calculated directly via the staggered spin-spin 
correlations at the T point. However, since the total magnetization Stot = SpA + Srs commutes 
with the Hamiltonian of the system, [8^*, H] = 0, the uniform spin gap A„ cannot be extracted 
from the uniform spin-spin correlations at the T point in a canonical quantum Monte Carlo sim- 
ulation. Instead, one obtains A„ = limk_s>r Au(k) from measurements performed at the finite 
momenta closest to the F point for each finite system. Supplementary Fig. [21 shows the finite size 
data for these gaps at U/t = 3.3, 4 and 4.5. For U/t = A, both gaps scale to finite values in the 
thermodynamic limit, with A^ being about four times as large as A^. For the other two values of 
U/t, both gaps clearly vanish in the thermodynamic limit. 
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4 Density correlations 



The density-density correlation function is given by 

where a,b E {A, B}. At half-filling, (rzxa) = ('^xh) = 1- Supplementary Fig. [3] shows the finite 
size scaling of the density correlations at the largest distance, C^'"^ = \C^'^{0, d^)! and C^'^ = 
|C^'^(0, di)| at U/t = 4. Both scale to zero in the thermodynamic limit, and no long-range 
density correlations persist. Furthermore, in comparison with the spin correlations, the density 
correlations are seen to be significantly weaker and essentially zero within the statistical error for 
system sizes L > 9. Consistently, we also find no long-range density ordering when analyzing the 
density structure factor (not shown). 

5 Dimer-dimer correlations - charge sector 

In this section, we present our results on the dimer-dimer correlations in the charge sector. We mea- 
sured both the correlations between the kinetic energy bond operators and the current operators. 
The spin sector is treated in the following section. 

Correlations between the kinetic energy bond operators 

/c(x, a; y, b) = ^{cl^^Cyf,^ + 4;,„c^aa), (10) 

a 

and the current operators 

j(x, a; y, b) = -i ^(^^^Cy^^ - ^^^c^a^) (11) 

a 

can be defined between both nearest-neighbour and next-nearest neighbour sites on the honeycomb 
lattice. 
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To probe for VBC order in the kinetic energy sector, we measured the three inequivalent 
dimer-dimer correlation functions 

Ci=\{{k{0,A-0,B)k{dL,A-dL + d„B)))\, 1 = 1,2,3, (12) 

at the largest distance on the finite lattices. For an illustration of the different relative bond 
orientations, see the inset of Suppl. Fig.lH The upper panel of Suppl. Fig. |4] shows the finite size 
scaling of the Cl at U / 1 = 4. These correlations scale to zero in the thermodynamic limit, hence 
no long-ranged bond order in the kinetic energy persists. Furthermore, in comparison with the spin 
correlations, these correlations are also seen to be significantly weaker. 

To probe for the persistence of nearest-neighbour currents in the ground state, we measured 
the current-current correlation functions between the bonds of the honeycomb lattice 

q = I ((j(0, A- 0, B)j{dL, A; dL + S„ B)))\, 1 = 1,2, 3, (13) 

at the largest distance on the finite lattices. The corresponding finite size scalings are shown in 
the lower panel of Suppl. Fig. HI Again, long range correlations in the thermodynamic limit can be 
clearly excluded, indicating the absence of currents between nearest neighbour sites in the ground 
state. 

To probe for bond order and currents between next-nearest neighbour sites, we measured 
all inequivalent next-nearest neighbour bond-bond and current-current correlation functions at the 
largest distances both within the same sublattice and between the two sublattices. 





= \{{k{S),A-v^,A)k{dL,A-dL 


+ r.,A)))|, 


(14) 




= \{{j{Q,A-T^,A)j{dL,A-dL 




(15) 


r^AB,i 


= \{{k{Q,A-vi,A)k{dL,B-dL 


+ r.,i?)))|, 


(16) 


^AB.i 


= \{{j{<d,A-v^,A)j{dL,B-dL 


+ v„B)))\, i = l,...,Q. 


(17) 
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For an illustration of the different relative bond orientations, see the inset of Suppl. Fig. [51 The 
quantum Monte Carlo data for the correlations within the same sublattice at f//t = 4 are shown in 
Suppl. Fig. [51 Both C^^'^ and C^^'^ all scale to zero in the thermodynamic limit. The correspond- 
ing correlations between the two sublattices similarly decay to zero in the thermodynamic limit 
(not shown). Thus no bond ordering nor currents persist between next-nearest neighbour sites in 
the ground state atU/t = A. 

6 Dimer-dimer correlations - spin sector 

In the spin sector, we measured the dimer-dimer correlation functions 

A,,H = (((S.-S,-^)(Sfc-S,-^))) (18) 

where ij and kl are each nearest neighbour sites on the honeycomb lattice. The quantum Monte 
Carlo results for these correlations are shown and discussed in the main text. 

We furthermore measured correlations between the spin-current operators 

j,(x, a; y, b) = -i ^{-IT {cl^aCyta - clba^^aa) (19) 

a 

as well as the spin-bond operators 

A;,(x, a; y, b) = ^(-l)"(cL«Cyb« + ^fe^Cxaa), (20) 

a 

for next-nearest neighbour sites. We measured these correlations between all inequivalent pairs of 
next-nearest neighbour sites both within the same sublattice and between the two sublattices at the 
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largest distance on the finite lattices, 

= \{{h{0,A;r,,A)hidL,A;dL + r,,A)))l (21) 

Cf^' = |((j.(0,A;ri,A)j,(di,A;d^ + r„A)))|, (22) 

= \{{kMA;r,,A)ks{dL,B;dL + r,,B)))l (23) 

Ctf = \{{js{0,A;r,,A)UdL,B;dL + r„B)))\, z = l,..,6. (24) 

For an illustration of the different relative bond orientations, see the inset of Suppl. Fig. [6l Sup- 
plementary Fig. [6] shows the finite size scaling of the correlations within the same sublattice at 
U/t = 4. They all decay to zero in the thermodynamic limit. The corresponding correlations 
between different sublattices show a similar behavior (not shown). Thus no spin-bond order nor 
spin-currents persist between next-nearest neighbour sites in the ground state atU/t = A. 

7 Flux quantization measurement for superconductivity 

In the flux quantization measurement, we thread a magnetic flux $, in units of the flux quantum 
through the centre of a torus on which the electronic system lies. From the functional form of 
the ground state energy with respect to the threaded flux, £'o($/$o), we can distinguish between 
normal and superconducting ground states. The signature of the latter requires that the macro- 
scopic energy difference Eo(^/^q) — _E'o($/$o = 1/2) scales in the thermodynamic limit to a 
periodic function of period 1/2, and the occurrence of an energy barrier between $/$o = and 
$/$o = 1 /2. In contrast, a metallic phase is characterized by Eq{^ /$o) - Eq{^ /<I>o = 1 /2) van- 
ishing as a power law as a function of system size, while in an insulating phase, it would vanish 
exponentially. Figure |7] compares the QMC results of the macroscopic energy difference at f/ = 
with that at t//t = 4. In both cases one clearly observes the vanishing of this quantity in the ther- 
modynamic limit. Hence, no signal for superconductivity is obtained from these flux quantization 
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measurements. 



8 Order parameters for superconductivity 

Order parameters for superconductivity are in principle obtained from considering the irreducible 
representations of the Dq point group of the honeycomb lattice, which can be described as a trian- 
gular lattice with a basis of two atoms in the sublattices A and B, respectively. The Cooper pair 
wave function of a superconducting state is a product of a spin, orbital and a sublattice component. 
Since Pauli's principle requires the wave function to be antisymmetric under particle exchange, 
we obtain the following possibilities for spin-singlet pairing: for an even (odd) orbital part, the 
wave function must be symmetric (antisymmetric) under sublattice exchange. It is convenient to 
introduce pair creation operators in the singlet channel 

At(x, a; y, b) = - i^c^^^, (25) 

where a, 6 G {A,B}. The operator 

At(x) = l(At(x, A; X, A) + At(x, B; x, B)) (26) 

describes on-site s-wave pairing, which is symmetric under sublattice exchange. In Suppl. Fig.[8l 
we show the s-wave pair-pair correlation function Ca = |(^s(0)Aj(dL))| at the largest distance 
at U/t = 4. No long-ranged pairing correlation sustains to the thermodynamic limit; instead, the 
on-site pair-pair correlation function decreases rapidly. 

Extended pair creation operators based on nearest neighbour pairing can be expressed in 
terms of phase factors /f , /2 , /s , « £ {^,B}, 

3 

At(x, ff, ft fi, /f , , fi) = J][/^At(x, A; X + 5„ B) + /f At(x, 5; X - (5„ A)] (27) 

i=l 
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For an extended s-wave, 

AL,„,(x)=At(x,l,l,l,l,l,l). (28) 
Nearest neighbour p-wave states relate to 

At^(x) = At(x,0, +1,-1,0,-1,+!), (29) 
At^(x) = At(x,0,+1,+1, 0,-1,-1), (30) 

and nearest neighbour ci-wave states to 

Alj^) = At(x,0, +1,-1,0, +1,-1), (31) 
Al^ ^(x) = At(x,-2,+l,+l,-2,+l,+l). (32) 

X —y 



In terms of next-nearest neighbours, one furthermore obtains the singlet /-wave state 

6 

A^(x) = ^(-1)' [At(x, A; X + r„ A) - At(x, 5; x + r„ B)]. (33) 



In order to probe for superconductivity based on nearest neighbour or next-nearest neighbour 
pairing in the above symmetry sectors, we directly measured in real-space the various inequivalent 
pair-pair correlation functions at the largest distances on the finite lattices, i.e. 

= |((A(0, A; 0, B)A\dL, A; + 6„B)))\, i = 1, 2, 3, (34) 

for the nearest neighbour pairing states, and 

C^^'^ = |((A(0,A;ri,A)At(di,A;di + r,,A)))| (35) 
C^^'* = |((A(0,A;ri,A)At(di,5;d^ + r„5)))|, 2 = 1,..., 6, (36) 
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for next-nearest neighbour pairing states both within the same sublattice and between the two 
sublattices. As shown in Suppl. Fig. |9] and Suppl. Fig. \Wi both nearest neighbour and next- 
nearest neighbour pair-pair correlation functions are very weak, even reaching zero within sta- 
tistical errors for L > 9. From this, we exclude pairing in all the above symmetry sectors, since 
the full Cooper pair correlations (Ae.,._.(0)AL.-.(di)), (A^MKMl)), (Ap.(O)At^(dz.)), 
(Arf,^(0)A;^^(di)), (Arf^,_^,(0)Ai^^_^^(dL)), and (A/(0)A}(di)) are linear superpositions of 
the above pair-pair correlation functions, and hence vanish in the thermodynamic limit. We can 
thus exclude superconductivity in the half-filled Hubbard model on the honeycomb lattice. 

9 Lower bound for singlet states in the RVB phase 

We consider the internal energy as a function of the projection parameter 9, 

^(0) ^ <^-ic"'-t;!"r'^^> ■ (37) 



where is the trial wave function (cf. Methods). Let {\n)} be the set of eigenstates of H. 
Then, we have 



E{e) -Eo = y e-(^"-^«)®(E„ - Eo)^P^'' 

^ ^ 1 + me) ^ ^ °M 0^ 



|2 



where Eq is the ground-state energy, |0) the ground-state, and define 
Next, we consider an energy scale 



(38) 



(39) 



e < (El - Eo)/N , (40) 

where Ei is the energy of the first excited singlet with the same quantum numbers as the ground- 
state, and define 6*, such that for 9 > 9*, (£'(9) - Eq)/N < e. Typically, e can be taken of the 
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order of the statistical error in the energy, such that the condition on e is clearly fulfilled. However, 
e is not defined in terms of the statistical errors; the only defining condition on e is (|40l ). With such 
a definition we have 

such that 



< 



n>0 



so that it also holds that R{Q*) <^ 1. The last inequality also implies that 



since the sum in i?(9) consists of positive definite terms. In case the overlaps in the last inequality 
are finite, 

{E, - Eo)e* » 1 , (44) 

such that 1/6* provides a lower bound for El -Eo- Incase | (l|^r) P/| (0|^r) P < 1 such that the 
inequality (l43l) is fulfilled due to a vanishing overlap, we miss the lowest excited singlet state, and 
1/6* provides a lower bound for the next lowest singlet with a finite overlap with the trial wave 
function. 

For the determination of the lower bound for singlet states we concentrated on the value 
U/t = 4, centered in the RVB phase. We verified that in the case L = 2, where the system 
can be fully diagonalized, setting e = lO^^t, 1/6* = (0.74 ± 0.04)t (s.e.m.) is a lower bound. 
Namely, for this system size, Ei — Eq = 1.84t for the first excited singlet state of momentum 
k = 0. We also verified that this state has an overlap of 0.22 with the trial wave function. The 
above value of e corresponds to the maximal error for all system sizes. The uncertainty in the 
determination of 6* is taken as the maximum between (i) the distance from 6* to the value of 6 
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for £"(6) /N = E{Q*) /N — e and (ii) the distance between two consecutive values of 6 around 6*. 
By means of error propagation, we then estimate the error in 1/0*. In Suppl. Fig.[TT]we display 
1/6* for L = 3, 6, 9, and 12, and an extrapolation to the TDL. In all these cases, the lower bound is 
above the spin-gap, as well as the extrapolation to the TDL. Hence, we find no evidence for singlet 
states with the same quantum numbers as the ground-state, that may become degenerate with it in 
the thermodynamic limit. However, this result does no exclude the possibility of low-lying singlet 
states that have vanishing overlaps with our trial wave function. 
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Supplementary Figure 1 | Comparison of the finite size scaling between the spin corre- 
lations C^'^ and Cf^^ at the largest available distance and the staggered structure factor 
Saf at U/t = 4 (upper panel) and U/t = 4.5 (lower panel), using 3rd order polynomials in 

1/L. Error bars denote standard errors. 
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Supplementary Figure 2 | Comparison of the finite size scaling between the staggered 
spin gap and the uniform spin gap A„ at U/t = 3.3,4 and 4.5 (top to bottom). The 
extrapolated values in the thermodynamic limit for U/t = A are = 0.023 ± 0.007 (s.e.m.) 
and A„ = 0.099 ± 0.001 (s.e.m.). Error bars denote standard errors. 
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Supplementary Figure 3 | Finite size scaling of the density correlation functions C, 



A,A 



and Cf^ a\U/t = 4. Error bars denote standard errors. 
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Supplementary Figure 4 | Finite size scaling of the nearest-neighbour correlations Q 
and Cj aXU/t = 4. The inset illustrates the three inequivalent directions with respect to 
the reference bond marl<ed by AB. Error bars denote standard errors. 
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Supplementary Figure 5 | Finite size scaling of the next-nearest neighbour correlations 
C^"^' and Cf^'^ at the largest distance at Vjt = 4. The inset illustrates the inequivalent 
directions with respect to the reference bond marked by AA, with the lower (upper) star 
corresponding to equal (different) sublattices. Error bars denote standard errors. 



38 



0.05 
0.04 
0.03 
0.02 
0.01 
O 

0.04 
0.03 
0.02 
0.01 

'^0 0.04""""0.rR5 Ti.l2 O.fS 0.2 0.24 0.28 0.32 0.36 

1 /L 

Supplementary Figure 6 | Finite size scaling of the next-nearest neighbour correlations 
and Cf:^f at the largest distance at U/t = 4. The inset illustrates the inequivalent 
directions with respect to the reference bond marked by AA, with the lower (upper) star 
corresponding to equal (different) sublattices. Error bars denote standard errors. 
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Supplementary Figure 7 | The energy difference of Eo{^/%) - E{^/^o = 1/2) for dif- 
ferent system sizes at f//t = and U/t = 4. Note that the scale for U/t = 4 \s four 
times smaller than for u/t = 0. The flattening of the energy differences exclude the su- 
perconducting ground state at both U/t = and U/t = 4. Error bars denote standard 
errors. 
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Supplementary Figure 8 | Finite size scaling of on-site s-wave pairing correlation Ca at 

U/t = 4. 
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Supplementary Figure 9 | Finite size scaling of nearest neigbour pair correlation at 
U/t = A. The inset illustrates the three inequivalent directions with respect to the reference 
bond marl^ed by AB. Error bars denote standard errors. 
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Supplementary Figure 10 | Finite size scaling of the next-nearest neighbour pair corre- 
lations C^^'^ and C^^'^ aXU/t = 4. The inset illustrates the inequivalent directions with 
respect to the reference bond mar[<ed by AA, with the lower (upper) star corresponding 
to equal (different) sublattices. Error bars denote standard errors. 



41 




1/L 

Supplementary Figure 11 | Lower bound 1/9* for L = 3, 6, 9, and 12 at f//t = 4. For 
comparison, the values of the spin-gap are reproduced. The inset shows the internal 
energy E{Q) as a function of the projection parameter 9 for the L = 12 system. The 
vertical bar gives the position of 9*. Error bars in and E{Q) denote standard errors. 
For the determination of the error bars in 1/9*, see the text. 
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